Abstract 



The calculation of a segment of eigenvalues and their corresponding eigenvectors of a Hermitian matrix 
or matrix pencil has many applications. A new approach to this calculation based on a density matrix has 
been proposed recently and a software package FEAST has been developed. The density-matrix approach 
allows feast's implementation to exploit a key strength of modern computer architectures, namely, 
multiple levels of parallelism. Consequently, the software package has been well received, especially in the 
electronic structure community. Nevertheless, theoretical analysis of FEAST has been lagging and that 
a convergence proof has yet to be established. This paper offers a detailed numerical analysis of FEAST. 
In particular, we show that the FEAST algorithm can be understood as the standard subspace iteration 
algorithm in conjunction with the Rayleigh-Ritz procedure. The novelty of FEAST is that it does not 
iterate directly with the original matrices, but instead iterates with an approximation to the spectral 
projector onto the eigenspace in question. Analysis of the numerical nature of this approximate spectral 
projector and the resulting subspaces generated in the FEAST algorithm establishes the algorithm's 
convergence. This paper shows that FEAST is resilient against rounding errors and establishes properties 
that can be leveraged to enhance the algorithm's robustness. The analysis here suggests that FEAST can 
be extended to non-Hermitian problems. Further investigations into numerical quadrature rule suitable 
for approximating spectral projectors are also worthwhile. 
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1 Introduction 

Solving matrix eigenvalue problems are crucial in many scientific and engineering applications. For problems 
of moderate size, robust solvers are well developed and widely available fTI and are sometimes referred to 
as direct solvers [?]• These solvers typically calculate the entire spectrum of the matrix or matrix pencil 
in question. In many applications, especially for those where the underlying linear systems are large and 
sparse, often only selected segments of the spectrum are of interest. A new approach for these calculations for 
Hermitian matrices and matrix pencils based on density matrices has been proposed recently Ql^ . From an 
algorithmic point of view, this new approach tries to compute an exact invariant subspace - which is the action 
of the density matrix - approximately, as opposed to, for example, the well-known Krylov subspace methods 
(see for example [H |6l [TH [20] ) which try to compute approximate invariant subspaces (the Krylov subspaces) 
exactly. The density-matrix approach maintains a basis for a fixed-dimension subspace but updates it per 
iteration. In this view, it is similar to the non-expanding subspace version of an eigensolver based on trace 
minimization |261 127j but with a different subspace update strategy. From an implementation point of view, 
this new approach is similar to spectral divide-and-conquer [31 [3] in that the calculation is expressed in terms 
of high-level building blocks that can much better exploit the advantages of modern computing architectures. 
In this case, the high-level building block is a numerical-quadrature based technique to approximate an exact 
spectral projector. This building block consists of solving independent linear systems, each for multiple right 
hand sides. A software package FEAST i22j based on this approach has been made available since 2009. 
Nevertheless, theoretical analysis of FEAST has been lagging its software development. In particular, a 
convergence proof of FEAST has yet to be established, and subtle numerical idiosyncrasies (that are spotted 
to arise occasionally [8]) lack explanations. 

This paper shows that the FEAST algorithm can be understood as a standard subspace iteration in con- 
junction with the Rayleigh-Ritz procedure, an approach which is explained in standard references such as 
[7], page 157, or [24], page 115. The important variation is that the subspace iteration in FEAST is carried 
out on an approximate spectral projector obtained by numerical quadrature. Our analysis shows that the 
quadrature approximation perturbs the projector's eigenvalues but not the eigenvectors. Consequently, the 
convergence of subspace iteration and hence of FEAST can be established by the same approaches shown 
in [24]. By exploring the structure of the generated subspaces, the numerical characteristics of FEAST can 
be much better understood so that all previously reported idiosyncrasies now have satisfactory explanations. 
Furthermore, FEAST's robustness can be enhanced by new techniques made possible by these structure anal- 
yses - techniques such as estimating the number of eigenvalues present or judging if the dimension chosen for 
the subspace in question is appropriate. This paper puts the FEAST algorithm on a more solid foundation, 
and we believe the overall analysis would allow us to generalize FEAST to non-Hcrmitian problems as well. 

*peter. tang@intel.com 
tpolizzi@ecs.umass.edu 
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2 Overview 



Consider two n-hy-n Hermitian matrices A and B, with B being positive definite: 

A:^A", and B^C"C, C invertible, 

where AI^ stands for complex-conjugate transposition of a matrix M. The generalized Hermitian eigenvalue 
problem (GHEP) is to determine eigenvalues and eigenvectors A and x, respectively, where 

Ax = XBx. 

The theoretical properties of GHEP are well known. In particular, the set of eigenvalues A are real valued, 
and one can choose a set of eigenvectors that are i?-orthonormal. In other words, let A be a diagonal matrix 
consisting of the n (real) eigenvalues, counting multiplicities. Then the eigenpair {X, A) satisfies 

AX^BXA, and X"BX = I. 

In particular, 

X-^ = X"B, and A = BX AX-^ = {BX)A{BX)". 

Given A and B and an interval on the real line I = [A_, A+], our goal is to obtain all of the eigenvalues that 
lie inside [A_,A+], and their associated eigenvectors. Our strategy is motivated by the spectral projector 
onto the invariant subspace corresponding to the eigenvalues of interest. This is the operator XxX^ B, 
where Xx is the set of eigenvectors spanning that invariant subspace. If we can compute {XxX^ B)y for any 
n-vector y, then applying XxX^ B on a large enough set of random vector^ Y = [yi, 2/2, • ■ • j J/p] will likely 
lead to 

spim{XxX^BY) = span(A:i), 

which is true as soon as rank( A^|^i?y) ~ rank(Xx). Once a basis Q for the invariant subspace is known, 
the target eigenpairs can be obtained, for example, by solving the (much smaller) p-hy-p GHEP of (j4p. Bp) 
where Ap — Q^AQ and Bp = BQ. This GHEP can be solved by standard solvers such as those available 
in LAPACK. 

Both the operator XxX^ B and the action [XxX^ B)y can be represented as a complex contour integral 
(more details later). Replacing the integral with a numerical quadrature rule yields an approximation formula 
for computing [XxX^ B)y for any y. To compute this numerical quadrature, it suffices to solve a number 
of linear systems with y as the right hand side. Furthermore, it turns out that span(Xi) is also an invariant 
subspace of the quadrature-based operator. Therefore, the straightforward subspace iteration algorithm 
applied with this quadrature-based operator can be used to capture the invariant subspace span(Ari). The 
eigenpairs of interest are then determined by a Rayleigh-Ritz procedure, which is the solution of a (much 
smaller) GHEP obtained by projecting A and B onto the captured subspace. The following is the flow of 
our analysis. 

1. We review the integral representation of XxX^ B and analyze the properties of the quadrature-based 
operator. In particular, we show that not only span(A'i) is an invariant subspace, but this subspace 
also corresponds to highly dominant eigenvalues. 

2. Based on the properties of the quadrature-based approximate spectral projector, we present a subspace 
iteration algorithm adapted for our original GHEP. We prove convergence of this iterative process in 
the sense that the distances between each of the eigenvectors x in Xx and the generated subspaces go 
to zero at a rate as fast as gap*"', where k is the iteration number and gap ^ 1 in practice. 

3. The structure of the subspaces generated in the iterative process are further analyzed to show that the 
approximate eigenvalues converge at a rate of gap^'^ while the corresponding residual vectors converge 
to zero at a rate of gap*". Furthermore, we show how the actual number of eigenpairs in [A_, A+] can 
be estimated robustly. 

4. Finally, we present a number of numerical experiments to illustrate the theoretical analysis. 

^ The effectiveness of randomized methods have been studied rigorously in a large body of recent works. Excellent surveys 
can be found in [l2lfT8] . 
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3 Computing Approximate Spectral Projection 



In this section, we consider the spectral projector XxX^ B as a function of matrix, representable as a contour 
integral^. We will show that when this integral is approximated by a numerical quadrature, the resulting 
approximate spectral projector has a number of useful properties. 

Let (X, A) be the eigenpair of the GHEP AX = BXA. Given an interval I — [A_,A+] on the real line, 
let C be the circle centered on the real line and intersecting it at exactly A_ and A_|_. Let 7r(A) be the 
complex- valued function defined by the contour integral (in the counter clockwise direction) 



ZTTL 



1 



A 



■ dz. 



\<ic. 



(1) 



The Cauchy integral theorem shows that 7r(A) = 1 for A inside the circle C and 7r(A) = for A outside of C. 
The idea is that applying tt to a matrix whose eigenvalues are A corresponds to a spectral projection. Since 
the GHEP is specified by two matrices, we define one single matrix that serves as a "surrogate." This is 



A = XkX-^ = XKX"B 

whose eigenpair is exactly (X, A) by design. Assuming just for now that neither A_ nor A+ is an eigenvalue 
of A, then the function 'n{A) is well defined as 

7r(i) = — i {zI-A)~^dz. 
2ni Jc 



On the other hand, 
Hence 



7r(i) = XTr{A)X-^ = Xtt{A)X"B = XiX§B. 
1 



XxX^ B 



2ttl 



{zl - A) dz 



as long as {A_, A+} does not intersect with A's spectrum. Moreover 

{zl - Af^ ^ [zI - XAX"Bf^ = [B-\zB - BXAX"B)f^ = {zB - A)'^B. 



(2) 



Hence the spectral projection of a set of n- vectors Y — [yi , ^2 , • ■ • , 2/p] admits an integral representation 
involving A and B: 



{XiXiB)Y (zI-A) Wdz^^ 

Zttl Jc Zttl 



(zB-A) ^BYdz. 



(3) 



Now that {XxX§B)Y is expressed as an integral, it is natural to approximate it via numerical quadraturesH. 
To this end, consider the scalar function 7r(A) in Equation [T] and the parametrization 4){t)^ —l<t<2>: 

=c + re'-5(i+*), and = e'^ , 

where c = + A_) and r = 5(A+ — A_). For A e R while not equal to either A+ or A_, 



^(A) = 



2ttl y_i (j){t) - A 



1 

2tu. 

1 

2tu 



1 - A 



dt, 



dt 



<P{t) - A (j){t) - A 



^4>{2-t)-X 
dt. 



dt 



^ This technique is well studied. The interested reader can find further details in 114) . 
^ See [25] for a different application of numerical quadrature to eigenvalue problems. 
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The integrand is a function of t as well as A. For a fixed A, a quadrature rule can be used to produce a value 
p(A) that approximates 7r(A). We now employ one fixed quadrature rule for all A G R\{A_, A+} and consider 
the results as a function p(A) of A. 

For a g-point quadrature rule on [— 1 . 1] (for example Gauss-Legendre) where the weights and nodes are 
(wfe, tfe), Wk > and — 1 < t/j < 1, fc = 1, 2, . . . , g, we have p{X) w 7r(A) where 

27rt^y 0(^fe)-A (t){tk)-Xj <t>k-X (j)k-\J 

4'k = (t>itk) and ^fc — u'fc(/)'(tfc)/(27rt). Furthermore, we will employ quadrature rules where |tfe| 7^ 1 for all k, 
which implies that p(A) is well defined for all A G K. As a result, the matrix function p{A) is well defined 
for all A, as its spectrum is real, including the case when either A_ or A+ is an eigenvalue. 
Applying Equation [3] to A, we have 

piA)Y = Y,U^kI-A) Y + Y,U^kI-A) Y. 

k=l k=l 

Applying Equation [5] gives 

p{A)Y = Yl=iU<t^kB-A)-'BY + YX=MT^B-A)~^BY, in general, 

= 2 X;Li {^^k{<f>kB - Ay'^BY^ , if A, B, and Y are real valued. ^ ' 

Equation [5] in general involves solving 2q systems of linear equations, each with p right hand sides. Only q 
systems need to be solved, however, if the GHEP is in fact real valued. Note however that since and (j)k 
are complex valued, the linear systems are complex valued. 

Intuitively, the function p{A) would approximate the spectral projector 7r(A). Let us analyze how well this 
approximation works. As a function of matrix, 

p{A) = p{XKX-^) = Xp{K)X-^ ^ Xp{K)X"B. (6) 

This shows that p{A) preserves A^s eigenvectors and hence the invariant subspaces while changing each 
eigenvalue from A to p(A). In particular, given any vector y expanded in the eigenvectors X: y = X^asa '^^^>'^ 

P(^)2/= I]«AP(A)XA. (7) 

AeA 

Equation [7] motivates our key strategy: try to obtain the invariant subspace of p{A) (instead of working with 
A directly) that corresponds to the eigenvalues p(A) for X E I. See Section [01 for an numerical illustration 
of Equation [71 

Turning our focus to the scalar function p(A), A € R, we note that it suffices to examine the reference case 
of [A_, A+] = [—1/2, 1/2] because p(A), due to our parametrization, satisfies 



Equation 3] shows that 

Prct(A) 



p(A) = p,ef((A-c)/(2r)). 

1 -A 1 - 2Acos(7r(l +tfe)/2) 

2 ^ l-4Acos(7r(l+tfe)/2)+4A2 ' 



Figure [T] shows the Picf(A) of an (7 = 8 Gauss-Legendre quadrature in normal as well as logarithmic scale. 
We see that components of y (cf. Equation[7]) associated with A G [—1/2, 1/2] will be attenuated by no more 
than a factor of 1/2. In fact, except for |A| very close to 1/2, those components remain basically invariant. 
On the other hand, components associated with |A| > 0.75 will be attenuated by at least a factor of 10"''. 
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Gaussian Quadrature with 8 Points ^ Gaussian Quadrature with 8 Points 




Figure 1: Gauss-Legendre quadrature using 8 points on [—1,1]. The left is PicfW in normal scale, and 
the right is |prcf(A)| in logarithmic scale. The picture on the left conveys the general idea that prcf(A) 
approximates 1 inside [A-,A+] = [—1/2,1/2], and 0, outside. The picture on the right illustrates how 
[prcf (A)| gets very small very soon beyond the reference [A_, A+] of [—1/2, 1/2]. 



Note that the range of Prcf (A) is essentially [0, 1] although the function can go "out of bound" by 0.02 or so 
above 1 or below 0. Throughout this paper, we will use the notation 1+ > jOiof(A) > 1/2 for A e [—1/2, 1/2] 
to say that Pref(A) is essentially in [1/2, 1] for A € [—1/2, 1/2]. The log-scale plot in Figure [1] shows that 
Prof (A) crosses zero a number of times for |A| > 1/2 and in particular Pref(A) is not monotonic in the region 
A > 1/2 or A < —1/2. Referring to Figure [TJ we can infer, for example, p{A)y will leave all components of 
y associated with A G [A_,A+] almost invariant but attenuate by at least four orders of magnitude those 
components associated with A farther than (A+ — A_)/4 from [A_,A+]. We will utilize Gauss-Legendre as 
our default quadrature rule. Prompted by Figure [U we were able to provcQ that for all Gauss-Legendre 
quadrature rules, /Orcf(A) > 1/2 for |A| < 1/2 and that piof(±l/2) = 1/2. 

4 Subspace Iteration on Approximate Spectral Projector 

The previous analysis shows that the approximate spectral projector is of the form p{A) — XTX^ B where 
F — /9(A). The invariant subspaces of p{A) are identical to those of A. Moreover, A's eigenvalues inside 
I — [A_, A+] are mapped to the dominant values in F, as 1+ > p{I) > 1/2 and |p(A)| < 1/2 for A ^ [A_, A+]. 
In fact, for A outside of [A_,A+], \p{X)\ <^ 1 except when A is quite close to [A_,A+]. We now employ a 
numbering convention to the eigenvalues and eigenvectors Aj, Xj using the ordering of the p{Xj)- We number 
the 7S such that 

71 > 72 > • • • > 7i > 1/2 > |7i+i| > • • • l7n|, 
and order the As accordingly so that 

7j=p(Aj), Axj^^XjBxj, p{A)xj^jjXj, j = l,2,...,n. 

In particular, i is the number of eigenvalues in I. These Ai, A2, . . . , Ai and their associated i?-orthonormal 
eigenvectors Xx = [xi, X2, ■ ■ ■ ^ xi] are our objects of pursuit. 

Subspace iteration is a standard pedagogical method (see for example [71111]) that can be used to capture 
invariant subspaces. While subspace iteration is seldom used in practice as the requirements for this method 
to succeed are stringent, the previous analysis on p{A) suggests it to be a perfect candidate for this simple 
iterative method to converge rapidly. This is because the invariant subspace of interest corresponds to highly 
dominant eigenvalues. The following is a standard form of subspace iteration but modified naturally to work 
for GHEP. In particular, orthonormalization is replaced with i3-orthonormalization. 

The proof is easy and elementary. 



6 



Algorithm SI (Subspace Iteration) 
1: Pick p random n- vectors Y(^q) = [j/i, 2/2, ■ • ■ i ?/p]- 
2: Set Q(o) ^ _B-orthonormalize(Y(o)). 

3: Comments: See explanation below on _B-orthonormalization. 
4: Set fc 1. 
5: repeat 

6: Yi^k) ^ P{A) ■ Q(k-l) 

7: Q{k) ^ i?-orthonornialize(F(fc)). 

8: fc -S- fc + 1 

9: until Appropriate stopping criteria 



B-orthonormalization is a generalization of orthonormalization. An n x p matrix Y is i?-ortlionormal iff 
Y^ BY = Ip. Consider now an n x p matrix Y of full rank. The followings are some elementary properties 
relevant to subsequent discussions. 

1. There exists a p x p upper triangular matrix R such that Q ~ Y R^^ is _B-orthonormal. This R can be 
obtained via the Cholesky factorization Y^ BY = R^ R. 

2. Another way to B-orthonormalize Y is to solve the p-hy-p GHEP 

AZ = BZA 

where 

A = Y"AY, B = Y"BY, 

for eigenvalues A and B-orthonormal eigenvectors (of length p) Z. As a result, Q = YZ is B- 
orthonornial. Z is not upper triangular in general, but this does not contradict the previous point. 

3. We say two _B-orthonormal matrices Q and Q' are equivalent iff there is a p x p unitary matrix U such 
that Q' = QU. Equivalent S-orthonormal matrices are B-orthonormal basis for the same subspace. 

4. The _B-norm of a vector is a generalization of the 2-norm: 

llyllij = Vy^By. 
If Q is an n X p _B-orthonormal matrix and w is a p-vector, then 

\\Qw\\b = ||w||2. 

5. Given a p-dimensional subspace in C" spanned by a B-orthonormal basis Q, the _B-orthonormal pro- 
jector onto this subspace is QQ^ B. Given any vector y e C", QQ^ By is the vector in span((5) that 
is closest to y in i?-norm: 

- QQ" B)y\\B = min \\z-y\\B. 

26span(Q) 

Some basic convergence properties of Algorithm SI are immediate consequence of well-established charac- 
teristics of subspace iterations (see for example [U [23l [24l HI]). Nonetheless, the fact that Algorithm SI 
utilizes an approximate spectral projector leads to properties not shared by general subspace iterationf[f|. 
First, application of p{A) to Y (Step 6 of Algorithm SI) often yields a set of almost-linearly-dependent vec- 
tors. Second, unlike a standard setting of subspace iteration, it is unrealistic to assume that the subspaces 
span((5(fc)) can capture all the p eigenvectors {xi, X2, ■ ■ ■ , Xp}. However, since the backbone of Algorithm SI 

^ Readers aequainted with acceleration methods can also view approximate spectral projection as a special acceleration, 
which warrant further analysis beyond the readily available results for plain subspace iteration alone. 
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is subspace iteration, our analysis below resembles that for this well-known method. We will highlight the 
differences as we proceed. 

Our analysis involves examining sections of the eigenvectors X = [xi, a;2, ■ . ■ , x„] and the corresponding 
eigenvalues. We set up some simplifying notations: For integer 1 < < n, 

Xg' — [xi+i,Xi+2, ■ ■ ■ ,Xn], 

Ti = diag(7i,72, . . . ,7^), 

r^' = diag(7i+i,7£+2,---,7n)- 

The main property of Algorithm SI is that the distance between each of the eigenvectors Xj, j = 1, 2, . . . ,p, 
and the generated subspace span((5(fc)) converges to zero at a rate |7p+i/7j|'^. This is favorable to p{A) as 
1+ > 7j > 1/2 for the eigenvectors of interest, j = 1, 2, . . . , i, while are mostly very small elsewhere (see 
Figured]). We describe these convergence results more formally in Lemma [T] and Theorem [T] below. These 
are straightforward generalizations from using 2-norm to that of i?-norm in Theorem 2 of |23) or Theorem 
5.2 on Page 118 of M. 

Lemma 1. Let Q be an n x p B-orthonormal matrix of the form 

|7p| > 0, where W is {n — p) x p and Z is p x p. Carry out one iteration in Algorithm SI: 
. Y ^ p{A)Q, 

• Q ^ B-orthonormalize(Y) . 
Then 

Q = {Xp + Xp,W) ■ Z-^ 

where 



ip+i 



'jh, j = l,2,...,p. 



7j 

Wj and Wj being the j-th columns of W and W , respectively. 
Proof. First note that 

p{A) = Xp Tp X^ B + Xp, Vp, X"^, B. 

Thus 

p{A)Q = {XpTp + Xp-Tp^W)- Z-\ 

= (Xp + Xp,w)rp ■ Z'\ w^rp,wr-\ 

Q = {Xp + Xp,W)Tp ■ Z-^ ■ R-\ for some R, 
= {Xp+Xp,W)-Z-\ 



Since the j-th columns of W and W are related by 

1 



Wj^—Tp,w.j, j = l,2,...,p, 

l3 



|7p+l| > |7p+2| > • • • > \jn\, 



as claimed □ 



7p+i 



Wwjh, j = l,2,...,p. 



The following establishes the convergence property of subspace iterations on p{A). 
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Theorem 1. Consider Algorithm SI. Suppose Y(o) such that the px p matrix BY(^q) is nonsingular. 
Then, 

1. there are constants aj,j — 1,2, ... ,p, and 

2. each Q^^-j is of the form 

Q^k)^{Xp + X^p,W)-Z-\ 

for some {n ^ p) x p matrix W and p x p matrix Z (the W and Z in general do not remain the same 
as iteration proceeds ), 

so that the j-th column, wj, ofW at the k-th iteration satisfies 

k 



\\wj\\2 < aj 



lp+1 



j = l,2,...,p. 



In particular, for each iteration k, and each j — 1,2, . . . ,p, there is Sj G span((3(fc)) such that 

k 



7p+i 



Proof. Since 



7i 



/ — XX^ B — XpX^ B + Xpi X^ B , 



Y(fi) — [XpXp B + Xp,X" B)Y(^a), 
= {Xp + Xp,W)X^BY^o), 



where 



W^iX^BY^o))iX^BY^o)) ■ 
After B -orthonormalization, (3(o) is of the form 

Q(o) = [Xp + Xp,W) ■ Z-\ 

Define 

"j = l|Wjl|2, j = l,2,...,p, 
as the 2-norm of the j-th column ofW. Applying Lemma\^ repeatedly shows that Qf^k) is of the form 

{Xp + Xp,W)- z-\ 

for some W where the j-th column, Wj, of W satisfies 

k 

aj, i = l,2,...,p. 



Ip+i 



7j 



Finally, let Sj be the j-th column of Q(^k^ Z. Thus Sj G span((5(j.)) and Sj = Xj + XpiWj, implying 



ll-'^j -XjWs = \\Xp,Wj\\B = \\wj\\2 < aj 



Ip+i 



This completes the proof □ 



Applications of the quadrature-based approximate spectral projection involve solutions of linear systems (see 
Equation [5]). We examine the implications when the solutions are inexact, producing quantities of the form 

p{A) Q + A, A is a n X ]3 matrix, 

whenever the projection is applied to a Q matrix in Step 6 of Algorithm SI. The error A accounts for the 
error produced in solving the systems 



(8) 
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(k) 

for Zj where k ranges over the quadrature points and j ranges over each column qj of Q. Each cohimn 
of A is the sum of forward errors over all quadrature points of the linear-system solutions for the column 
in question. We assume these forward errors are of order roundoff 0{u) with the constant dependent on 
the condition numbers of C (where B = C^C) and that of 0^/ — A (see [13] chapters 7-9 for example). 
Lemma[2]and Theorem [5] below are the modifications to Lemma[T]and Theorem [T] that take A into account. 
They show that Algorithm SI is resilient to errors in the linear solver used to implement p{A). 

Lemma 2. Let Q be an n x p B-orthonormal matrix and m be an integer m < p. Suppose Q is of the form 

Q^[{Xra+Xp,W + Xra'F)-Z~^ V]-U 

where W is (n — p) x m, F is (n — m) x m, Z is m x m, V is n x (p — m), and U is unitary of dimension 
p X p: U^U = Ip. Carry out one iteration in Algorithm SI where the linear systems are solved inexactly: 

• y -s— p{A)Q + A (assume Y remains full rank p), 

m Q -(r- B-orthonormalize{Y) . 

Define the partition 



X-^AU 



H 



Z 
/ 



[ -^m I -^rn' ] — 



(top) 
m 

(bot) 



E„ 



Em is n X m, Em' is n x (p~ m), and i?,n°^'' is the rax ra top portion of Em- and is the (n — m) x m 

bottom portion of Em ■ Suppose r,„ + E^n"^^ is nonsingular, then Q is of the form 

iXm+Xp,W + Xm'F)-Z-^ V]-U, 
W , F, V , Z, U are of same dimensions as W, F, V, Z, U , respectively, and U is unitary. Furthermore, 

w = Tp,wirm + El^°p'^)'\ 

F = Tm'FiTm + + E^^^'HTm + Et"^)''. 

Denote (Tm + Em°^"') by (/ + . Key relationships between W, F and W, F are 





< 


7p+i 






7i 




< 


Ip+i 






7m 


\F\\2 


< 


7m+l 






7m 



IICIhllM^IU, J = 1,2,, 



(1 + IIC||2)|1W^||2, 

(l+llCI|2)||F||2 + ^(l + ||CI|2)||i?(r'|l. 
l7m| 



We make a few comments before presenting the proofs. A typical choice of p has the property of |7p+i| ^ 1. 
A natural by-product is that |7j| will also be small for some of the js, j < p, that are close to p. The m in 
the next Lemma, while formulated as a general integer m < p, is meant to correspond to the index for \^m\ 
when \^m\ is not too small. This m exists because we know 1^ > 71 > • • • > 7i > 1/2. For such m, the Z 
and Z matrices that arise during Algorithm SI remain well conditioned (in fact these matrices tend to the 
identity matrix) as will be evident later on. Consequently, the magnitude of Em will remain tame, that is, 
its norm will be not much bigger than that of A. 

Proof. Similar to the proof of Lemma\^ 

p{A)Q = [{XmTm + XpiTpiW + Xm'^m' F)-Z-' p(A)V ]-U, 

A - [{XmE^^°P^+Xm'E^^°'y)-Z-^ XEm' ] ■ U, 

Y = [{Xm{rm+Ei^°P'>)+Xp,rp,W + Xm'{Tm'F + E^:°'^))- Z-^ V ]-U, 

= [{Xm+Xp,W + Xm'F)-{Tm+Et''^)-Z-^ V ]-U, 
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where 

w = rp,VF(r,„ + s(^°p))"\ 

B-orthonormalize X„i + XpiW + X„iiF so that {X„i + XpiW + Xm' F)Z^^ is B-orthonormal. Complete a 
B-orthonormal basis in span(y) to obtain a B-orthonormal basis of the form 

{Xm+Xp,W + X^,F)Z~^ V 

Hence Q must be equivalent to the above, meaning that there is a p x p unitary matrix U such that 

(X^ + Xp,W + Xm'F)Z~^ VyU. 

w = rp,w(i + o^;n\ 
F r„-F(/ + c)r,-i + £;i^°')(/ + c)r,-i. 



Finally, let (/ + T:^^E, 



-1 p(top)N 



Consequently, 



Ip+i 



Ip+i 



7p+i 



WChWh, j = l,2,...,m, 



7m 

7p+i 



7r) 



(1 + IICI|2)||W-||2, 

(l + ||Cl|2)||F||2 + ^||4r^||2(l+||CI|2). 
l7m| 



WWjh < 
\\W\\2 < 
\\Fh < 

This completes the proof □ 

Theorem 2. Consider Algorithm SI. Suppose Y(o) such that the p x p matrix X^ BY^Q'j is nonsingular. 
Let the application of the quadrature-based projector be of the form 

p{A)Q + A 

as in Lemma \^ Suppose there is an integer m, m < p and a constant k such that for all iterations the 
matrices Em°^^ and Em"^"^ defined as in Lemma [H satisfy the followings: 

\\{I + T-^E(l°^'>y' -Ih<n, and |l4'°*^||2 < ac. 
Then there are constants a.j,j = 0,1,2, ... ,m, and that Q(^k) can be represented as 

Qik) = [ Xrn + Xp,W + Xrn'F V ] U, 

U being p x p unitary. Moreover, for each iteration k and each j = 1,2, ... ,m. 



II/.II2 < 



k(1 + k) 

l7m| 



fe-1 



dcf 



7j 



In particular, for each iteration k and each j = 1,2, . . . ,m, there is Sj £ spaii((5(fc)) such that 



fe-i 



\sj - XjWb < aj Pp^j + ao /3^;J ^(1 + n) 



£=0 



f ^ k{1 + k) 

l7m| 



^((l + «;)/3„,„f. 



i=0 
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Proof. Since 



where 

Hence Yj^o-) is of the form 



I — XX^ B — XpXp B + Xpi X^, B, 



= {X.pXp B + Xp,X" B)Yi^Q), 
= {Xp + Xp,M)X^BY^o), 

M=iX^BY^o)){X^BY^o))~'- 



Y^o)^[X^ + Xp,W + X,n'F V]-G 

where W is the first m columns of M and F is simply the zero matrix of dimension (n — to) x to. Now 
B-orthonormalize Xm + XpiW + Xm'F into {X^ + Xp/W + X^'F) ■ and complete a B-orthonormal 
basis for span(Y(o)) to yield 

[{Xra+Xp,W + X^^,F)Z-^ V]. 

Hence Q(q) must be equal to 

[iXrn + Xp,W + X^,F)Z-^ V]-U 
for some p x p unitary matrix U . Define the constants aj, j — 0,1,2, ... ,m, 

"0 = \\W\\2, aj = \\wj\\2, 

as the 2-norms of the matrix W as well as each of the individual columns. Apply Lemma\^ repeatedly shows 
that Q(^k) of the form 

Q(k) = [X„,+ Xp,W + Xm'F V ]-u. 

The W and F (as well as V and U, though we do not care so much about these two) matrices change as 
iteration proceeds while maintaining the relationships 

WWjh < l3p,j\\Wj\\2+ l3p,rnK\\W\\2, 
\\W\\2 < Pp,m{l + ^)\\W\\2, 

k{1 + k) 



\\F\\2 < I3rn^rnil + H)\\F\\2 + 



Simple recurrence analysis show that ||VF||2 o,nd \\F\\2 (note that in the beginning \\F\\2 = 0) at the k-th 
iteration satisfy 

\\W\\2 < aoPlmil + K)\ 

\\F\\2 < ^^j±{^J2{{l + n)P^,mY. 

17m I 

With the bound just obtained on |1W^|12 cit each iteration, the simplification jSpj < (3p^m, and simple recurrence 
analysis, we obtain bound of ||wj||2 o,t the k-th iteration as 

fe-i 



\\wj\\2 < ajpl^ + ao/3^J ^ (1 + k)^ 



£=0 

Finally, since 

= [ Xm + Xp,W + Xra'F V ] ■ U, 

the j-th column of Xm + Xp/W + Xm'F is in span((5(fe)). This j-th column is 

= Xj + Xp' Wj + Xm' fj 
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and 



\\Sj ~Xj\\b 



— 1 1 Wj + X„i> /j 1 1 _B , 

-^m' f j\\B ^ 

= l|wj||2 + ll/jib, 



1=0 



k(1 + k) 

\lm\ 



This completes the proof □ 

Let us compare Theorems [T] and by considering the eigenvectors of interest, Xj, j — l,2,...,i. TheoremU] 
guarantees that as long as p > i, the eigenvectors wih be approximated as close as one wishes by the 
generated subspace provided one iterates long enough. In practice, however, one would aim at p large 
enough so that |7p+i/7i| ^ 1 to ensure fast convergence. Note that fast convergence, or mere convergence, 
can be expected only for those j , 1 < j < m where |7p+i/7m| is reasonably small. It is therefore possible 
that the p-dimensional subspaces only carry m < p "dimensions" of useful information. In the situation 
where |7m+i| ~ l7m+2| « • • • « |7p| « l7p+i| *C 1, the directions oi p — m dimensions of Q(k) are more or 
less a random artifact of rounding errors; although is well-conditioned (to the extent of the condition 
of C) due to being i?-orthonormal. It is thus crucial to show that the "first" m > i dimensions of Q(k} will 
evolve robustly regardless of the potential non-convergence of the remaining p ^ m dimensions. Starting 
with Theorem [21 our analysis assumes only that |7p+i/7m| <C 1 for some m < p and allow implicitly that 
the remaining p — m dimensions, usually denoted by the letter V, to be arbitrary. Compared to Theorem [1] 
Theorem [2] also produces upper bounds on the distance between Xj and the generated subspaces. In contrast, 
these upper bounds do not go to zero. First, the term aj/S^j is the same as that in Theorem [TJ Next, as 
long as K is moderate, 0{u) or even 0{\/u), u being machine roundoff, the second term goes to zero as 
|7p+i/7m| ^ 1 for some m > i. Finally, as long as (1 -I- K)/3m,m < 1, the third term remains bounded. This 
means that error in linear system solution does not affect convergence in any fundamental way. It limits 
the eventual accuracy to a threshold commensurate with the linear solver's accuracy, of which k serves as a 
quantifier. Even if (1 -I- K)(3m.m > 1, it is close to 1 as k is small in general. The last term is approximately 
kK/l'yml- As long as both f3pj and /3p,m <C 1, a few iterations will allow the eigenvectors of interest be 
approximated to the level of accuracy of the solver before kK/\jm\ grows large. 

In summary, subspace iteration has a self-correcting nature. Errors introduced in one iteration are attenuated 
to some extent in subsequent iterations. Errors are somewhat limited only to those incurred in the latest 
iteration. Section 16.21 presents numerical illustrations of Theorem [2] 

5 Eigenproblems 

The previous section shows that if the subspace dimension p in Algorithm SI is chosen large enough so 
that |7p+i/7i| <C 1, then the generated subspaces — span(Y(fc')) — span((3(j,')) will capture rapidly the 
eigenvectors xi,X2, . . . ,Xi. In fact, if |7p+i/7m| ^ 1 for some m, i < m < p, the subspace will also capture 
the additional eigenvectors a:i+i, Xi+2, ■ • • , a;,„ very well. This scenario is typical when there are eigenvalues 
outside of [A_, A+] but close to the boundaries. This means that 7i « 7i+i « • • • « 7^ for some m > i. In 
order for |7p_|_i/7i| <C 1, a p > m has to be chosen. As a by-product, more eigenvectors then the i targeted 
ones will be captured. The analysis below considers Algorithm SI using a subspace dimension p presumed 
to be big enough so that |7p+i/7m| ^ 1 for some integer to, i < m < p. 

Now to complete the story, we must show how to make use of these subspaces that have presumably captured 
the wanted eigenvectors, to actually obtaining them. Specifically, we wish to obtain the xi,X2, ■ ■ ■ ,Xm and 
the associated Ai, A2, . . . , A„i. Intuitively, these to eigenvalues will be among the p eigenvalues of the reduced 
problem Aqz = Xz, Aq = Q"AQ or that of the GHEP Ayw = XByw, Ay ^Y" AY and By =Y" BY. 
The corresponding eigenvectors will be among the p vectors Qz or Yw, respectively. We now validate this 
intuition. 

The first theorem examines the structure of the subspaces more closely. The desired properties about the 
eigenpairs of those reduced problems will follow naturally. 
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Theorem 3. Let Q be a n x p B-orthonormal matrix of the form 

Q = [ {Xrn + X^,Wra') ■ V] 

for some m < p where the 2-norm of each column of Wm' , Wi is less than e for some < e <C 1 

\\wi\\2 < e, i = 1, 2, . . . , m. 

Then the exist an equivalent B-orthonormal matrix Q to Q such that the n x p matrix BQ has the 
structur^ 



x"bq 



G 



G 



Al2 

A21 

Aoff + Adiag- 



An 



i22 



I„i is the m X m identity, G is {n — m) x (p — to) orthonormal, that is, G^G — Ip-m, and the As are 
small perturbations. The off-diagonal perturbation is of order e while the diagonal perturbation is of order 
€^ : ||Aoft||2 = 0{e) and ||Adiag||2 = 0{t^). The G, A12 and A22 terms are simply nonexistent for the case 
m = p. 

We comment here that the techniques used in estabhshing Theorems [3] and S] are very similar to those in |28) . 
In essence, our theorems can be considered as generahzing some of the results there from m — p to m < p. 

Proof. The last p — m columns V in the matrix Q can be expressed as 

V = XraW^m + Xr,,>Um' 

for some matrices Wm and Um' ■ Hence 

Q = [ {Xm + Xrn'Wm') ■ XmWm + X^'U^' ] ■ 

Observe that the m columns Xm -\- Xm'Wm> are close to being B-orthonormal because 

\\{Xm + Xra'Wm'f B (X„, + X^'Wrn') " Imh = WW^^W^'h = O(e'). 



Hence they can be made exactly B-orthonormal by a Z ^ chosen to be close to Im-' Z ^ = I„ 
IIA11II2 = O(e^), and {Xm + Xm'Wm')Z~^ being B-orthonormal. Next, 



A 



11: 



{XmWra + Xra'Um'f B [X^ + X^'Wra') Z' ' = 

implies + U^,Wm' = 0. Thus WWmh = 0{e). This fact, together with 

{XrnWm + Xm'Um') B [XmWm + Xm'Um') — Ip-m 

imply that \\U^,Um' ~ Ip-m\\2 — O(e^). This means that Um can be orthonormalized into G, Um = G R where 
\\R - Ip-mh = 0(6^). Hence Um' is of the form Um' = G + A22 where G"G = Ip-m and || A22II2 = O(e^). 



Putting all these together, 

x"bq 



B 



yH 

^m' 

Z-^ Wm 
Wm'Z-^ Um' 

G 
G 



{Xm+Xm'Wm')Z ^ XmWm + Xm'U„ 



A12 

A21 
Aoff + Adiag, 



All 



A. 



with ||Aoff||2 = 0(e) and ||Adiag||2 = O(e^) □ 



The displayed equation is not "drawn to scale." G in general is very tall and skinny. 
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Consider now the eigenvalues A of the original GHEP and partition A into 



A 



A, 



A, 



where A™ = diag(Ai, A2, . . . , Am) and Km' = diag(Am+i, Am+2, • • • , A„). The next theorem studies the p 
eigenpairs of the matrix AS where is a matrix of the form in the previous theorem. 

Theorem 4. Let A — diag(Am, Am') and S be an n x p matrix of the form 







G 


+ 



^12 



A21 



All 



\22 



U 



ioff 



^diag 



where Im is the m x m identity matrix, G is [n — m) x {p — m) orthonormal, G^G = Ip-m- That is, S 
is U perturbed by small off-diagonal and diagonal structures. Suppose ||Aoff||2 < e/(4||A||2) and ||Adiag||2 < 
e^/(4|| A||2). Then, therearem eigenpairs, (Ai, ii), (A2, 5;2), ■ • • , (Am, ^m), (among thep eigenpairs) of S^ AS 
such that ^ 

|A,-Aj|<e2 + J\ , j = l,2,...,m, (9) 

where rj is the gap between the spectrum of Am and G^A„i'G: rj = min|A — ii\ where A ranges over Xj, 
1 < j < m and ji ranges over all p ~ m eigenvalues of G^ A„i'G. And 

\\A{Sx,) ~ \,{S~x,)h < 2(e + e')(l + \\Ah/v). 
The values of rj in Equations\^ and 1 1 0\ are considered +00 in the case m = p 
Proof. It is straightforward to see that 

S"AS = U"AU 

= C/^AC/ + Soff + Sdiag 

where 



(10) 






E ' 




■ El 










+ 





E2 



Eos = Afff A?7 + C/^AAoff + A^^AAoff + Afff AA 



AfagAf/ + C/^AAdiag + A^^gAAdiag + AfffAAoff . 



By assumptions on the sizes of Aog and Adiag, we have H-Eofflb < e and ||-Ediag||2 < f^- 

Am 



U"AU = 



G^ Ajn'G 



whose p eigenvalues consist of the m eigenvalues of Am and thep — m eigenvalues of G^ Am'G, each of which 
lies in the range of Am' (because G^ G = Ip^m)- 

It is known that Hermitian perturbation of the form 

U"AU + E^if = 



Am 







E ' 


G Am' G 


+ 


E" 






admits a possibly tighter eigenvalue perturbation bound (see for example |P1 \19[ \29f } than that of Weyl's 
perturbation theorem (see for example [Tj', page 198). In particular, ive can use the elegant formula in 1171 
and conclude that there are m among the p eigenvalues A' of AU + Eos that are very close to the target 
eigenvalues: 

^ ' - + ^^7^ + A\\E,n\\'i ~ V + V^F+I^ 

where rj is the spectral gap between Am and G^ A„i'G. That is, rj min|A — /i| where A ranges over 
the spectrum of Am and ji ranges over the spectrum of G^ Am'G. In essence, when rj is not too small. 
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\\j ~ I < O(e^), instead of the bound of \Xj ^ A^ j < 0(e) obtainable by simple application of the Weyl 
perturbation theorem. We note that the remaining p — m eigenvalues of AU + Eqs are also very close to 
the p — m eigenvalues A^'G, but this fact does not concern us too much. 

To complete the bound on the eigenvalues, we apply Weyl's theorem to the i?diag perturbation to AU + Eog 
and conclude that there are m eigenvalues Xj , j ^ 1,2, . . . ,m, of AS such that 

\X,-X,\<e^ + , „ j = l,2,...,m. 

ri + \/ri^ + 'it 

It is obvious that \Xj — Xj \ < should m — p. Or for convenience, we adopt Equation\^with rj = +00 should 
m — p. 

Moving on to the eigenvectors, we use (A,i) to denote any one of the m eigenpairs {Xj,Xj) of S^AS. First, 

\\U"AUx - Xi\\2 < \\E,ff + £;diag||2 < e + e'. 
Partition x into the upper m elements a and lower p — m elements b, we have 

\\U"AUi- Xxh > \\{G"A„,,G - A/)6||2 > vPh, 

implying \\b\\2 < {e + e^)/r]. 

\\AUx~XUxh < |lA„,a- Aa|l2 + |1A„,,(G&)- A(G6)|l2, 

< (e + e2) + 2||A||2||&||2, 

< {e + e') + 2{e + e')\\Ah/v. 

Finally, 

\\ASx-XSS;\\2 < \\AUS:-XUx\\2 + \\AAx-XAS:\\2, 

< 2(e + e2)+2(e + e2)||A||2/7y, 

< 2{e + e'){l + \\A\\2/v)- 

It is obvious that b is nonexistent should m — p. In other words, rj can be taken as +00 in Equation \10\ 
.should m — p. This completes the proof □ 

Theorems [3] and m together show that subspace iteration apphed to the approximate spectral projector p{A) 
can be used to obtain target eigenpairs of the original GHEP Ax = XBx. A high-level reasoning it that 
Algorithm SI with an appropriate dimension p > i will generate subspace basis Q with property satisfying 
Theorem [3l Theorem |4] then says that target eigenpairs can be well approximated by {X,Qx) where (A,x) 
are among the eigenpairs of the reduced system AQ. The next theorem provides the details. 

Theorem 5. Let Algorithm SI be carried out to a certain iteration so that Q satisfies the requirements in The- 
orem[3[for an integer m > i. Then the p eigenpairs of AQ can be numbered (Ai,yi), (A2,y2), • • ■ , (Am,2/m) 
so that 

\Xj-Xj\ and \\A{Qyj)- X,B{Qyj)\\2, j = l,2,...,m 
are small as stated in Theorem^ 

Proof. By Theorem\^ Q is equivalent to another B-orthonormal basis Q such that BQ is of the form 

G 



X"BQ = S = 



A = U + A, 



G^G — Ip-m- Because Q and Q are equivalent B-orthonormal basis, there is apx p unitary matrix V such 
that Q = QV, V"V = Ip. 

Q"AQ = V"Q"AQV ^ V"{S"AS)V. 
By Theorem^ there are m eigenpairs {Xj, Xj), j — 1,2, . . . ,m of S^ AS such that 

\Xj - Xj\ and \\A{Sxj) - Aj(5x~ )||2 
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are small. Clearly, {Xj,yj) = (Xj Xj) are the corresponding eigenpairs ofQ^AQ. We already know that 
the XjS approximate the target eigenvalues. We now show that Qijj have small residuals. 

\\A{Qy,)-~X,B{Qy,)\\2 = \\A{Qi,) - X,B{Qi,)\\2, 

= \\BXA{X"BQ)ij~X,B{Qi^)\\2, 

= \\BXASxj~XjB{Qij)\\2, 

= \\X-"x"BXASS:j - XjX-" X" B{Qij)\\2, 

= \\X-"{KSxj-XjSxj)\\2, 

< \\C\\2\\ASxj-XjSij\\2, 

where B = C^C and X being B-orthonormal must be of the form X — C^^W for some unitary matrix W . 
This completes the proof □ 

One final detail before we formally state the FEAST Algorithm. Q is the result of B-orthonormalization of 
Y in Algorithm SI. As stated previously, one way to B-orthonormalize Y is to solve the (dimension p) GHEP 
defined by the pair {Y^ AY, Y^ BY), which is equivalent to solving AQ had we known Q. A most natural 
way to augment Algorithm SI is therefore to solve the GHEP (Y^AY, Y^ BY) at every iteration. This serves 
two purposes at once: to B-orthonormalize Y , and to obtain approximate eigenpairs of the original GHEP 
[A, B) as iteration proceeds. This is the essence of the FEAST Algorithm [21] stated below in the language 
of this current paper. 



Algorithm FEAST 
1: Pick p random n- vectors ^(0) ~ [j/i , 2/2 , ■ • • , J/p] • 
2: Set (5(0) ^ i3-orthonormalize(Y(o)). 
3: Set ^ 1. 
4: repeat 

5: Y(fc) ^ p{A) ■ Q(k-1) 

6: Form A ^ Y^^ AY^^) and B ^ F^f^ BYt^k) ■ 

7: Solve the p-hy-p generalized Hermitian eigenvalue problem (GHEP): AX — BXA. 

8: g(fc) ^ Y(fe) • X. 

9: Comments: Q{k) is -B-orthonormal. Part of A and Q(k) approximate the target eigenpairs. 
10: k ^ k + 1 

11: until Appropriate stopping criteria 



A number of comments are in order. 

1. Suppose p is chosen so large such that some of the |7j|s for j < p are very small. It is obvious then that 
Y(i) would be numerically rank deficient. The direct solver employed in Step 6 will likely complain as 
the Cholesky factorization of B will fail. The straightforward way to prevent this failure is to perform 
a rank-revealing QR factorization [TO] and replace Y(i) with a smaller set of fewer columns (a smaller 
p). In practice, the Cholesky factorization in LAPACK [T would reveal the row index at which failure 
occurs. The FEAST implementation in |22| employs the simple resizing technique of using only the 
columns of Yi before failure. This convenient, though less robust than rank-revealing QR, method 
works well in practice. 

2. Suppose p > i is chosen and in fact there is an m, i < m < p such that \jp+i/jm\ ^ 1 while |7p+i/7p| 
is not too small. Theorem |4] shows that as iteration proceeds, m of the p eigenvalues obtained at each 
iteration will approximate the actual eigenvalues of the original GHEP very well. In particular, we will 
very likely see i of these fall inside [A_, A+] and the other m — i, outside. But as Theorem |4] suggests, the 
other p — m computed eigenvalues can in general be any average of the remaining n — m eigenvalues of 
the original problem. Indeed, some of these may very well fall inside [A_, A+J. In general, one does not 
know the value of i a priori. Consequently, one may not be able to judge unequivocally which among 
the p eigenvalues of the reduced system correspond to the actual targets inside [A_ , A+] . Identifying 
the i target eigenvalues from the p numbers can be made easier if i is known. An ideal situation is that 
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i is indeed known and that one observes exactly i of the p eigenvalues of the reduced system lie inside 
[A_, A-f.]. We will show in a moment how the value of i can be estimated in practice. 

3. If the subspace iteration is carried out with a p value that is actually smaller than i, then Algorithm 
SI will fail to converge almost certainly: Any good quadrature rule for computing the approximate 
spectral projection will map all interested eigenvalues to [1/2, 1] (or very nearly so); and in fact most 
of them will be mapped closely to 1. Algorithm SI will not converge in general as the ratio |7p+i/7i| 
will be very close to 1. 

4. Based on the previous discussions, the tasks of (1) estimating i in the case of p > i, and (2) judging if 
p is sufhciently large are important. We now discuss a method to accomplish both. 

The key is that the eigenvalues of the B matrices in the inner loop of Algorithm SI offer a robust estimate 
for the number i. Theorem |6] captures this point. 

Theorem 6. Let Q be the B-orthonormal basis in Algorithm FEAST at some iteration with 

Q=[iXm+Xra'W)Z-^ V], 

where each column ofW, Wi, is small: ||u'i||2 < e ^ 1. Consider Steps 5 and 6 of Algorithm FEAST: 

Y ^ p(i)Q, 
B ^ Y"AY. 

Then, up to 0{e) (and in fact most likely 0{e^)), i of B 's eigenvalues are inside [1/ A, 1^] and the other p — i 
are inside [0, 1/4). 

Proof. Following the analysis in Theorems \^ and^ Q — QU where 

X"BQ = S= ^ 

with U^U ~ Ip, G^G = Ip-m, A is of 0{e^) in the diagonal blocks and, 0{e), off diagonal. 

Y = p{A)QU, 

B = U" Q" {p{A)f Bp{A)QU, 

= U" Q" {XTX"Bf B{XTX"B)QU, 
= {Q" BX)TX"BXT{X"BQ)U, 

= u" s" s u. 



So the eigenvalues of B are the same as that of 



r2 , 

m' 



where G^G — Ip-m o,nd E is of 0{e^) on diagonal blocks and 0(e) off diagonal. (See Theoremsl3\and^for 
more details.) Considering E as a perturbation, the unperturbed matrix has m eigenvalues exactly equal to 

1+ > 7? > • • ■ > 7f > 1/4 > 7?+i >--->ll- 

Furthermore, because 7^ > 7^+1 > • ' ' ^ In, the remaining p ~ m eigenvalues are in the range [0,7^ 
The differences between the perturbed and unperturbed eigenvalues are bounded by an expression of the form 
+ 2e^/(77 + yjrf + 4e2) where rj is the spectral gap between the r^„ and G^ F^, G. This completes the proof 

□ 



m+l] 
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Based on Theorem [6l computing i?'s eigenvalues offers a practical estimate of i. If p happens to be chosen 
appropriately in the first place, we find in practice that _B's eigenvalue count as early as fc = 2 already 
estimates i correctly. When B has no eigenvalues much smaller than 1/4, this is an indication that p needs 
to be adjusted larger. Eigenvalues or singular values of other matrices involving p(j4), i?, Y, or Q can also 
be used to estimate Tm and hence offer a way to estimate i. We discuss two other instances at the end of this 
section. Nevertheless, using i?'s eigenvalues is natural as the matrix is readily available. Algorithm FEAST 
with Estimate incorporates this estimation step. 

Algorithm FEAST with Estimate 

1: Pick p random n- vectors Y(o) — [j/i, 2/2, ■ • • , J/p]- 

2: Set Q(o) -i— _B-orthonormalize(y(o)). 

3: Set k ^ 1. 

4: repeat 

5: ^ p{A) ■ Q(k-l) 

6: Form A ^ Y^^^ AY^^) and B 4- Y^^^ BY^^) ■ 

7: Solve the p-hy-p generalized Hermitian eigenvalue problem AX — BXA. 

8: Q(k) ^ Y(^k) ■ X. 

9: Comments: Q^k) is -B-orthonormal. Part of A and Q^k) approximate the target eigenpairs. 
10: if k equals 2 then 

11: Compute B's eigenvalues; count how many are > 1/4; also examine the smallest one. 

12: Estimate i. If p is inappropriate, adjust p, goto Line 1 and start over. 

13: end if 

14: Comments: If fc > 2, then an estimate of i is available. Check the i smallest residuals corresponding 

to computed eigenvalues in [A_, A+]. 
15: fc ^ fc + 1 

16: until Appropriate stopping criteria. 



For completeness, there are other matrices whose eigenvalues or singular values can serve as an estimate of i, 
the number of targeted eigenvalues inside [A_, A+]. Consider a i?-orthonormal matrix Q that approximates 
the target eigenvectors to a certain extent. Q ~ QU and X^BQ has the form of Theorem [3l Consider the 
matrix Z: 

Y ^ p{A)Q, 

z ^ q"by. 

Then the number of Z's eigenvalues no less than 1/2 would be an estimate of i. The reason is 

Z = U"Q"BXrx"BQU, 

= u"{s"rs)u, 



= U 



H 



E U. 



Finally, consider the case when the original problem is a simple eigenvalue problem, HEP. That is, B = I. 
Let Q be an orthonormal matrix that approximates the target eigenvectors to a certain extent. Consider the 
matrix Y where Y <— p{A)Q. Then the number of Y^s singular values no less than 1/2 would be an estimate 
of i. The reason is 



Y 



XTX^QU, 



= XT 



X 



G 



Fm' G 



U + E, 
G]U + E. 



This is a perturbation to [XmTm XmiTmiG] U whose singular values are that of [XmTm XmiVm'G] (as U 
is unitary), which are in the range of |F| and exactly i of which are in [1/2, 1"*"]. 
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6 Numerical Experiments 



We present a number of numerical experiments to illustrate various aspects of our previous analyses. To this 
end, we utilize synthetic, controlled, examples as opposed to eigenproblems that arise from actual applica- 
tions. Given the problem dimension n and a search interval [A_,A+], we generate A, the diagonal matrix 
containing the eigenvalues, somewhat randomly, except for special placements of some of the eigenvalues 
near the boundaries of [A_,A+]. Random unitary matrices are the basic ingredient of our test matrices. 
With a specified condition number k, random matrix C is generated as UY,V^ where U and V are random 
unitary matrices and E are random singular values so as to make the condition number of C equal k. The 
matrix B is constructed as i? = C^C. The eigenvectors X are constructed by solving CX = W where W 
is a random unitary matrix. Finally, the matrix A is constructed a.s A = {BX) A (BX)^ . This way, 

AX = BXA, 

and {X, A) is the eigenpair of the generalized eigenproblem defined by A and B. 
6.1 Approximate spectral projector via quadrature 

A crucial property of the quadrature-based approximate spectral projector p{A) is that it preserves A's 
eigenvectors and changes only its eigenvalues from A to F = p{A) (see Equations [5] and [7]) : 

K 

p{A) ^a-(7/cS - A)-' ■ B = Xp{A)X"B = Xp{A) X'K 
fc=i 

We generate matrices A, B, A, X (as outlined previously) of dimension n = 300 with the elements of A to 
be uniformly distributed in [—30, 30]. The C matrices used in generating B — C^C have condition number 
100. We used Gauss-Legendre quadrature rule with 6, 8 and 10 quadrature points on [—1, 1]. For each test 
system and quadrature rule, we compute 

e =^ max e-j = p{A)xj - p{\j)xj . 

1<J<" ||A||2 

For each quadrature rule, 200 test cases are generated and Table [T] tabulates the maximum, mean, and 
standard deviation of these 200 es. 



Statistics of Quadrature Points of Gauss-Legendre 



{\\p{A)x,~p{\,)x,\\/\\A\\} 


6 




8 




10 




Maximum 


5.5 X 10- 


15 


9.0 X 10- 


15 


1.0 X 10- 


14 


Mean 


2.1 X 10" 


16 


2.4 X 10" 


16 


2.9 X 10- 


16 


Standard Deviation 


5.5 X 10- 


16 


8.0 X 10- 


16 


9.8 X 10- 


16 



Table 1: Key Property of Quadrature- Based Approximate Spectral Projector. For each eigenpair {\j,Xj) of 
A, we check if indeed p{A)xj p(Xj)xj. 



6.2 Convergence of subspace iteration using approximate spectral projector 

To illustrate Theorem [U we generate a complex generalized problem of dimension n — 500. We use Gauss- 
Legendre quadrature with 8 points on [—1, 1]. [A_, A+] is set to [15, 17]. The n eigenvalues are generated as 
follows. We pick four eigenvalues in [15, 17] by picking three randomly with uniform distribution in the region 
[15.2, 16.8]. The fourth is set to be 17. This guarantees that |p(Aj)| sa 1 for j = 1,2,3, and |p(A4)| = 1/2. 
These 4 eigenvalues are the only ones in [15, 17] and hence i = 4. Next, five eigenvalues are set to be in the 
interval (17,18] such that the values of |p(Aj)| are for £ = 3,5,7,9,11. The remaining 491 eigenvalues 
are chosen randomly with uniform distribution on the set [—40, 14] U [18,60]. The iteration of Algorithm 
Subspace Iteration is carried out with p = 8. With this choice of p, |7p+i/7j| is 2~^^ for j = 1, 2, 3, and 2"-'^°, 



20 



2~^, up to for the next 5 eigenvalues. Since the problem is generated, the eigenvectors Xj arc known, 
and the projectors P(fe) — Q{k)Qfk)^ are easy to compute. We examine the quantities ||(/ — P(k))xj\\B for 
each oi j = 1, 2, . . . , 8 for 5 iterations k = 1,2,. ..,5. Indeed these norms decrease in a way consistent with 
what the theorem predicts, except when the ultimate threshold of machine precision is reached. Table [5] 
tabulates the result. 









log2 1 




)Xj\\B 


at Iteration k = 


j 


logs 


79 




1 


2 


3 


4 


5 


1 


-11 


-12.9 


-23.8 


-34.7 


-43.3 


-43.2 


2 


-11 


-12.0 


-22.8 


-33.8 


-43.3 


-43.3 


3 


-11 


-11.6 


-22.4 


-33.4 


-43.4 


-43.4 


4 


-10 


-12.7 


-22.6 


-32.5 


-41.2 


-41.1 


5 






-7.3 


-15.1 


-23.1 


-31.1 


-38.5 


6 


-6 


-5.0 


-10.8 


-16.8 


-22.8 


-28.8 


7 


-4 


-3.2 


-7.0 


-11.0 


-15.0 


-19.8 


8 


r 




-1.1 


-3.0 


-5.0 


-6.9 


-8.9 



Table 2: Subspace Convergence, Complex CHEP. The convergence rate is consistent with the factor |7p-|_i/7j |. 
This test problem is designed with p — % and 7p+i = 2~^'^ . There are 4 eigenvalues in [A_,A+] = [15,17] 
and the nature of p{A) is that 1+ > 7^ > 1/2 for all 7^ e [A_,A+]. Convergence rate for the targets are 
hence quite uniform, unlike standard subspace iteration applied with the original matrices. 

To illustrate Theorem [21 we repeat the same experimental setting except that we added artificial errors 
to the linear solvers. To every solution z of equation of the form of Equation |S1 

^k{(i)kB - A)z = Bq, 

we modify z by a random error of 2^^u, u being the machine precision, 

z z + 2^^u \\z\\2 S, each element of the n-vector S is uniformly random in [—1/2, 1/2]. 

According to the bound of Theorem [21 which we restate here (see that section for details) 

lis, - x,\\b < a, P";.^ + ao /3p,™ ^(1 + k)^ + ^ + n)p^^„,)\ 

1=0 l"^™! i=o 

we expect the overall convergence rate not to be affected. Only the ultimate accuracy limit is degraded 
commensurate with the artificial errors injected here. The parameter m is flexible. Thus for each eigenvalues 
Xj, we can apply the bound with m = j. The bound suggests that the actual convergence limit is affected 
by the last term with the factor l/|7m|. Table [3] is consistent with these predictions. 

6.3 Eigenvalue and residual norm convergence 

We illustrate important aspects of Algorithm FEAST as stated in Theorem [5l The first example is complex 
GHEP, dimension 500, with [A_ , A-|_] — [15, 17]. We generate i = 5 eigenvalues well inside this interval. Eigen- 
values outside of [15, 17] are generated randomly except for a few specially placed so that 7 = 2~^-~^'~'^'~^ . 
Had p be set to 5 = i, the convergence rate would be somewhat slow. With p set to 8, convergence rate 
for the target eigenpairs will be linear with a factor of 2~^. The implication is that the three "collaterals" 
pair will also converge, except at a slower rate. This example reflects a typical scenario according to our 
experience with actual applications. There are often eigenvalues outside but quite close to the boundaries 
of [A_,A+]. As a result, the successful p will be strictly bigger than i and that the iterations will also 
obtain extra eigenpairs that can be called "collaterals." Table [H shows the numerical details. The ratios are 
|7p+i/7j| = 2~^ for the target eigenpairs. Note that eigenvalues accuracies improve by 2~^^ per iteration 
as suggested by Theorem [5l This is typical, especially when the collaterals converge. In this event, unless 
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j 


l0g2 


79 




2 


lo 

3 


ffn lie/- 

62 II V-* 

4 


5 


B at Iteration k 
6 7 


8 


9 


1 


-11 


-21.35 


-32.34 


-34.39 


-34.50 


-34.40 


-34.46 


-34.36 


-34.34 


2 


-11 


-22.97 


-33.66 


-34.34 


-34.38 


-34.33 


-34.37 


-34.33 


-34.38 


3 


-11 


-23.04 


-33.68 


-34.37 


-34.35 


-34.38 


-34.42 


-34.44 


-34.29 


4 


-10 


-19.51 


-29.51 


-33.40 


-33.35 


-33.39 


-33.44 


-33.39 


-33.40 


5 


-S 




-16.22 


-24.22 


-31.17 


-31.38 


-31.31 


-31.39 


-31.36 


-31.32 


6 


-6 


-12.20 


-18.21 


-24.21 


-29.16 


-29.31 


-29.43 


-29.45 


-29.34 


7 


-4 


-7.48 


-11.48 


-15.49 


-19.49 


-23.49 


-26.96 


-27.36 


-27.36 


8 


r 




-5.02 


-7.02 


-9.01 


-11.01 


-13.01 


-15.00 


-17.00 


-19.00 



Table 3: Subspace Convergence with Error in Linear System Solutions. Note that the overall convergence 
rate is not affected by errors injected into the solutions of linear systems. The ultimate accuracy achieved is 
consistent with the error bound of Theorem [2j By Iteration 8, the generated subspaces have captured the 
best they ever can the eigenvectors 1 to 7. The ultimate achievable accuracy degrades by a factor of 2 from 
eigenvectors 3 to 7, consistent with the factor of l/|7j|, j = 3, . . . , 7. 



Convergence of Eigenvalues and Residual 



log2 



log; 



|A,--Aj| 



at Iteration k = 
1 2 3 



-33.74 -51.15 -64.04 

-32.07 -49.50 -62.32 

-34.40 -51.88 -61.90 

-36.59 -54.23 -62.11 

-34.79 -52.18 -61.23 



\\Axj~\jBxj\\2 

at Iteration k — 



log; 



-28.57 -37.56 

-27.83 -36.82 

-29.23 -38.22 

-30.46 -39.46 

-29.58 -38.57 



-46.57 
-45.83 
-47.22 
-48.46 
-47.58 



-52.93 -52.95 

-52.95 -52.86 

-53.15 -53.10 

-53.14 -53.14 

-52.88 -52.91 



above are i target eigenvalues; below are "collaterals" 



6 


-6 


-30.73 


-40.64 


-52.61 


-24.92 


-30.91 


-36.91 


-42.91 


-48.54 


7 


-4 


-24.27 


-30.73 


-38.74 


-20.03 


-24.04 


-28.05 


-32.06 


-36.07 


8 


-2 


-24.55 


-28.95 


-32.99 


-19.96 


-22.08 


-24.08 


-26.09 


-28.10 



Table 4: Convergence of Eigenvalues and Residual Vectors. This table represents a typical scenario. Subspace 
dimension p is bigger than i but the "extra" dimensions also capture additional invariant subspaces, albeit 
slower. Note the eigenvalues converge linearly at the rate of (79/7^)^, while residuals do so at that of |79/7j|. 



the gap between the target and collateral eigenvalues are small, Theorem [5] predicts linear convergence of 
eigenvalues with the factor (jp+i/jj)'^. 

The second example is similar to the first: complex GHEP, dimension 500. We generate i = 5 eigenvalues 
well inside this interval. Eigenvalues outside of [15, 17] are generated randomly except for five specially- 
placed ones. One is placed so that 7 — 2'^, and four others are placed so that 7 is strictly bigger than, but 
extremely close to, 2~^. The other 491 eigenvalues are random but at least 0.5 away from [15, 17]. By setting 
p = 9, the convergence rate of the targets eigenvalues are expected to be linear with a factor (2~^)^ = 2~^^ 
or smaller. But the collaterals do not converge. Table [5] exhibits this phenomenon. 

The relationship between the subspace dimension p and the actual number of targets i can be subtle. In a 
typical scenario, p > i and that the collaterals will also converge, except at a slower speed. But in the case 
when the collaterals do not converge, one might think that there is no fundamental harm in carrying them 
along except for a moderate increase of computational cost. The analysis that accompanies Theorems 2] 
and [5] suggests some potential problems. Consider the previous example where the 9-dimensional subspaces 
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Convergence of Eigenvalues 



j 


log2 


79 

7j 




1 


lo^ 

2 


. |Aj-A,| 

ll^lb 
3 


at Iteration 
4 5 


6 


1 


-9 


-38.77 


-56.77 


-62.87 


-63.00 


-64.45 


-66.45 


2 


-9 


-35.38 


-53.39 


-61.17 


-62.87 


-62.65 


-62.55 


3 


-9 


-37.50 


-55.51 


-62.41 


-62.37 


-63.45 


-61.81 


4 


-9 


-36.77 


-54.78 


-65.55 


-65.45 


-63.13 


-63.00 


5 


-9 


-43.19 


-61.21 


-63.17 


-62.21 


-62.13 


-64.87 



above are i target eigenvalues; below are "collaterals" 



6 


-0.0023 


-33.53 


-35.55 


-35.56 


-35.56 


-35.56 


-35.56 


7 


-0.0017 


-32.36 


-37.61 


-37.62 


-37.62 


-37.63 


-37.63 


8 


-0.0012 


-31.26 


-36.75 


-36.77 


-36.77 


-36.77 


-36.77 


9 


-0.0006 


-30.29 


-35.07 


-35.07 


-35.07 


-35.07 


-35.07 



Table 5: Non- Convergence of Collaterals. There are 5 targets, and subspace dimension p is set to 9. The 
ratios |7p+i/7j| ss 1 for j = 6,7,8,9 and thus the collateral eigenvalues do not converge. These iterations 
would have been successful even if p was set to be just 5. 



capture the i target eigenvectors well, but not much of anything else. The reduced systems carry with 
them two subsystems. One is approximately Aj, and the other of the form Ai'G (in the notations of our 
theorems). If one is unlucky to have the eigenvalues of the second subsystem closely approximating some of 
the targets, convergence speed of target eigenvalues may be reduced to improvement of |7p+i/7i| per step, 
as opposed to the square of that. More importantly, some of the eigenvectors may actually be wrong! The 
residual may not converge to zero. We illustrate this phenomenon in the next example. For simplicity, we use 
a real- valued simple eigenvalue problem of dimension 500. We place just one eigenvalue A = 16 in the middle 
of [A_, A+] = [15, 17] but place two eigenvalues at 15-C and 17-f C so that p(15-C) = p{l7 + C) = 2"^. The 
remaining 497 eigenvalues are randomly generated except at least at a distance 3 away from [15, 17]. We set 
p = 2 and thus the target eigenvalue should converge at least by 2~^ per iteration, but usually at 2~^® per 
step. We contrivedly start the iterations with two vectors, one close to the target eigenvector, and the other 
about the middle of the two eigenvectors associated with 15 — C and 17 + (. That is, the Raleigh quotient 
with this vector is exactly 16. Table [6] illustrates the problem with a small gap between and G^A,n'G. 
As exhibited there, one of the two eigenvalues of the reduced system converge to 16, albeit only improving 
by 2^^ per step. Neither residual vector converges in any practical sense. 



p = 2, 



73 

71 

Sk is 



= 2- 



Convergence Hampered by Spurious Eigenvalues 
Examine log2((5fe/|| AII2) at Iterations k — 
2 3 4 5 6 7 



=1,2 |Aj - 16| 



nimj=i_2 



-18.45 -27.46 -37.60 -47.42 -46.00 -46.42 -46.00 -46.19 



-6.13 



-6.13 



-6.20 -12.66 -15.84 -16.13 -15.92 -16.46 



Table 6: 0{e^) Convergence of Eigenvalues and Non-Convergence of Residual. This artificial example is set 
up so that there is only one eigenvalue, A — 16, in the target interval. With p — 2 the ratio 7^-1.1/71 — 2~^. 
In fact, 72/71 = 2^^ as well. The collateral space is affecting the overall convergence. Convergence of 
eigenvalue falls back to 2~^ per iteration, not at the often enjoyed speed of 2~^® per iteration. More 
importantly, the residual vector is not really converging. The l/rj factor in Theorem]?] is realistic. For this 
example, convergence will be restored to the perfect situation had p be set to 1. 
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6.4 Estimation of eigenvalue count 

The number of eigenvalues in [A_ , A_|_] is valuable information; but guessing what that number is by counting 
the number of computed eigenvalues of reduced systems that fall inside [A_,A_|_] is an unsound practice. 
Theorem ini suggests that we can instead count the number of B's eigenvalues > 1/4. In the following 
complex GHEP example of dimension 48, we generated 8 eigenvalues inside [A_, A+] = [15, 17]. We place 
on each side of [A_,A+] 20 random eigenvalues of similar distribution to increase the chance of "spurious" 
eigenvalues. We set p to 12. Table [7] shows that the distribution of i?'s eigenvalues is a much more robust 
indication of i than that of computed eigenvalues of reduced systems. 



p ~ 12 eigenvalues, /ij, of B at Iteration k = 



3 


2 


3 


4 


5 


1 




1.031865 


1.042407 


1.042588 


1.042589 


2 




1.021825 


1.037967 


1.041980 


1.042530 


3 




0.999644 


1.000591 


1.000626 


1.000662 


4 




0.998992 


0.999999 


1.000000 


1.000000 


5 




0.998206 


0.999997 


1.000000 


1.000000 


6 




0.996776 


0.999936 


0.999999 


1.000000 


7 




0.929957 


0.999593 


0.999857 


0.999909 


8 




0.872044 


0.989169 


0.998930 


0.999845 


9 




0.201241 


0.210605 


0.211077 


0.211304 


10 




0.137805 


0.146882 


0.150190 


0.153307 


11 




0.086650 


0.095591 


0.098854 


0.104075 


12 




0.050975 


0.078311 


0.084910 


0.088754 


*H > 


1/4 


8 


8 


8 


8 


#A, e [A 


-,A+] 


10 


9 


9 


9 



Table 7: Eigenvalue Count of B to Estimate i. In this example, there are exactly i = 8 eigenvalues in 
[A_, A+] = [15, 17], p = 12. The computed eigenvalues of reduced problem may have more than 8 falling 
inside [A_, A+]. But the eigenvalue count of B estimates i correctly from Iteration 2 onwards. 



Along the same line, the next example in Table |5] shows that we can get an early indication that p is set 
too small by the eigenvalues of B. The example's setting is similar to the previous one, except p is set to 6, 
which is 2 less than the number of eigenvalues inside [A_, A+] = [15, 17]. The actual computed eigenvalues 
do not converge, which was to be expected. 





Computed ei^ 


;envalues Aj 






Eigenvalues /ij of B , 






of reduced system 


at Iteration k ~ 


as a monitor, at Iteration k = 


3 


2 


3 


4 


5 


2 


3 


4 


5 


1 


15.30888 


15.30960 


15.30709 


15.30358 


1.03473 


1.03600 


1.03768 


1.03997 


2 


15.53633 


15.54323 


15.54300 


15.54067 


1.01624 


1.01788 


1.02042 


1.02356 


3 


16.18678 


16.21928 


16.22838 


16.23129 


1.00012 


1.00016 


1.00016 


1.00016 


4 


16.54569 


16.58401 


16.58713 


16.58817 


0.99987 


1.00002 


1.00002 


1.00002 


5 


16.59844 


16.63769 


16.66165 


16.66953 


0.99871 


0.99984 


0.99984 


0.99985 


6 


16.81077 


16.83859 


16.85640 


16.86408 


0.74160 


0.89480 


0.96709 


0.99060 



Table 8: Eigenvalue Count of B to Judge p. In this example, there are i = 8 eigenvalues in [A_, A+] = [15, 17]. 
But p is set too small at p = 6. Computed eigenvalues will not converge in general. This table illustrates 
that a too-small-p can be detected by examining B^s eigenvalues as early as at the second iteration. The 
symptom is that none of S's eigenvalues are less than 1/4. 
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7 Conclusions 



We have shown that quadrature-based approximate spectral projectors are superb tools to be used with 
the standard subspace iteration method. This combination is the essence of the recently proposed FEAST 
algorithm and software ([HI [52]). Our detailed analysis establishes FEAST's convergence properties and 
show how its robustness can be further enhanced as methods for counting target eigenvalues and detecting 
inappropriate subspace dimension are identified. Eigenproblems of large- and-sparse systems fit FEAST 
naturally as it can tolerate less-accurate solutions of linear systems, allowing the use of iterative linear 
solvers (see Example 3 in [21]). The underlying analysis offers a natural approach for future work on the 
non-symmetric problems - bi- iteration ([S], or |30| page 609) with approximate spectral projection. The 
main tasks will be applying the quadrature-based projector in a form with left and right eigenvectors. While 
generalization of Theorems[3]and|3]will be needed, encouraging experimental results are already available [T5] . 
On a different note, we have used Gauss-Legendre quadrature as our numerical integrator of choice for /5(A) 's 
accurate approximation to the characteristic function 7r(A) on [A_,A4-] (see Equation [T]) . Nevertheless, 
accurate approximation of 7r(A) is by no means the only relevant property of an integrator suitable for 
FEAST. Investigation of other quadrature rules are worthwhile. One observation is that |p(A)| needs not 
approximate 1 very well on a large portion of [A_, A+] or decay to zero outside of [A_, A+] remarkably, both 
phenomena of which Gauss-Legendre possesses. It suffices to have, for example, p(A) fluctuates as long as 
1^ > p(A) > 77 ^ on [A_, A+] while keeping |p(A)| uniformly small outside [A_ — 6,X+ + 6] for some small 
6 > 0. Another observation is that it is valuable to have a quadrature rule that provides increasing accuracy 
by progressively adding more nodes (while maintaining the existing ones) . This would require us to find an 
alternative to Gauss-Legendre. In short, opportunities for further work abound. 
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